R Notebook

library(tidyverse)

── Attaching core tidyverse packages ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.4     ✔ readr     2.1.5
✔ forcats   1.0.0     ✔ stringr   1.5.1
✔ ggplot2   3.5.1     ✔ tibble    3.2.1
✔ lubridate 1.9.4     ✔ tidyr     1.3.1
✔ purrr     1.0.2     ── Conflicts ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the ]8;;http://conflicted.r-lib.org/conflicted package]8;; to force all conflicts to become errors

library(Seurat)

Loading required package: SeuratObject
Loading required package: sp
‘SeuratObject’ was built under R 4.4.0 but the current version is 4.4.1; it is recomended that you reinstall ‘SeuratObject’ as the ABI for R may have changed
‘SeuratObject’ was built with package ‘Matrix’ 1.7.0 but the current version is 1.7.1; it is recomended that you reinstall ‘SeuratObject’ as the ABI for ‘Matrix’ may have changed

Attaching package: ‘SeuratObject’

The following objects are masked from ‘package:base’:

    intersect, t

library(patchwork)
library(SingleR)

Loading required package: SummarizedExperiment
Loading required package: MatrixGenerics
Loading required package: matrixStats

Attaching package: ‘matrixStats’

The following object is masked from ‘package:dplyr’:

    count


Attaching package: ‘MatrixGenerics’

The following objects are masked from ‘package:matrixStats’:

    colAlls, colAnyNAs, colAnys, colAvgsPerRowSet, colCollapse, colCounts, colCummaxs, colCummins, colCumprods, colCumsums, colDiffs, colIQRDiffs, colIQRs, colLogSumExps,
    colMadDiffs, colMads, colMaxs, colMeans2, colMedians, colMins, colOrderStats, colProds, colQuantiles, colRanges, colRanks, colSdDiffs, colSds, colSums2, colTabulates,
    colVarDiffs, colVars, colWeightedMads, colWeightedMeans, colWeightedMedians, colWeightedSds, colWeightedVars, rowAlls, rowAnyNAs, rowAnys, rowAvgsPerColSet, rowCollapse,
    rowCounts, rowCummaxs, rowCummins, rowCumprods, rowCumsums, rowDiffs, rowIQRDiffs, rowIQRs, rowLogSumExps, rowMadDiffs, rowMads, rowMaxs, rowMeans2, rowMedians, rowMins,
    rowOrderStats, rowProds, rowQuantiles, rowRanges, rowRanks, rowSdDiffs, rowSds, rowSums2, rowTabulates, rowVarDiffs, rowVars, rowWeightedMads, rowWeightedMeans,
    rowWeightedMedians, rowWeightedSds, rowWeightedVars

Loading required package: GenomicRanges
Loading required package: stats4
Loading required package: BiocGenerics

Attaching package: ‘BiocGenerics’

The following object is masked from ‘package:SeuratObject’:

    intersect

The following objects are masked from ‘package:lubridate’:

    intersect, setdiff, union

The following objects are masked from ‘package:dplyr’:

    combine, intersect, setdiff, union

The following objects are masked from ‘package:stats’:

    IQR, mad, sd, var, xtabs

The following objects are masked from ‘package:base’:

    anyDuplicated, aperm, append, as.data.frame, basename, cbind, colnames, dirname, do.call, duplicated, eval, evalq, Filter, Find, get, grep, grepl, intersect, is.unsorted,
    lapply, Map, mapply, match, mget, order, paste, pmax, pmax.int, pmin, pmin.int, Position, rank, rbind, Reduce, rownames, sapply, setdiff, table, tapply, union, unique, unsplit,
    which.max, which.min

Loading required package: S4Vectors

Attaching package: ‘S4Vectors’

The following objects are masked from ‘package:lubridate’:

    second, second<-

The following objects are masked from ‘package:dplyr’:

    first, rename

The following object is masked from ‘package:tidyr’:

    expand

The following object is masked from ‘package:utils’:

    findMatches

The following objects are masked from ‘package:base’:

    expand.grid, I, unname

Loading required package: IRanges

Attaching package: ‘IRanges’

The following object is masked from ‘package:sp’:

    %over%

The following object is masked from ‘package:lubridate’:

    %within%

The following objects are masked from ‘package:dplyr’:

    collapse, desc, slice

The following object is masked from ‘package:purrr’:

    reduce

Loading required package: GenomeInfoDb
Loading required package: Biobase
Welcome to Bioconductor

    Vignettes contain introductory material; view with 'browseVignettes()'. To cite Bioconductor, see 'citation("Biobase")', and for packages 'citation("pkgname")'.


Attaching package: ‘Biobase’

The following object is masked from ‘package:MatrixGenerics’:

    rowMedians

The following objects are masked from ‘package:matrixStats’:

    anyMissing, rowMedians


Attaching package: ‘SummarizedExperiment’

The following object is masked from ‘package:Seurat’:

    Assays

The following object is masked from ‘package:SeuratObject’:

    Assays

library(celldex)

Attaching package: ‘celldex’

The following objects are masked from ‘package:SingleR’:

    BlueprintEncodeData, DatabaseImmuneCellExpressionData, HumanPrimaryCellAtlasData, ImmGenData, MonacoImmuneData, MouseRNAseqData, NovershternHematopoieticData

library(pheatmap)

# Speeds up scTransform
#BiocManager::install("glmGamPoi")

# Load AR and Tol datasets
AR_data <- CreateSeuratObject(Read10X(data.dir = "~/Desktop/R_projects/Dr_Krummey_Lab/heartTransplant_data/GSE262851_RAW/AR/"), min.cells = 3, min.features = 200)
AR_data$condition <-  "AR"
AR_data <- PercentageFeatureSet(AR_data, pattern = "^mt-", col.name = "percent.mt")

VlnPlot(AR_data, features = c("nCount_RNA", "nFeature_RNA", "percent.mt"))

Warning: Default search for "data" layer in "RNA" assay yielded no results; utilizing "counts" layer instead.

# I filtered a bit looser, nFeatureRNA > 5000 vs 2500, percent.mt < 15 vs 5
AR_data <- subset(AR_data, subset = nFeature_RNA > 200 & nFeature_RNA < 5000 & percent.mt < 15)
VlnPlot(AR_data, features = c("nCount_RNA", "nFeature_RNA", "percent.mt"))

Warning: Default search for "data" layer in "RNA" assay yielded no results; utilizing "counts" layer instead.

Tol_data <- CreateSeuratObject(Read10X(data.dir = "~/Desktop/R_projects/Dr_Krummey_Lab/heartTransplant_data/GSE262851_RAW/Tol/"), min.cells = 3, min.features = 200)
Tol_data$condition <- "Tol"
Tol_data <- PercentageFeatureSet(Tol_data, pattern = "^mt-", col.name = "percent.mt")
VlnPlot(Tol_data, features = c("nCount_RNA", "nFeature_RNA", "percent.mt"))

Warning: Default search for "data" layer in "RNA" assay yielded no results; utilizing "counts" layer instead.

Tol_data <- subset(Tol_data, subset = nFeature_RNA > 200 & nFeature_RNA < 5000 & percent.mt < 15)
VlnPlot(Tol_data, features = c("nCount_RNA", "nFeature_RNA", "percent.mt"))

Warning: Default search for "data" layer in "RNA" assay yielded no results; utilizing "counts" layer instead.

# Lets merge asap so we know the same transformations are applied
# Project is just a label to group the merged datasets and I never reference it. Feel free to change!
data.merge <- merge(AR_data, Tol_data, add.cell.ids = c("AR", "Tol"), project = "Heart")

# SCTransform instead of normalize -> variable features -> scale. More up-to-date
data.merge <- SCTransform(data.merge, vars.to.regress = "percent.mt", verbose = FALSE)
data.merge <- RunPCA(data.merge, verbose = FALSE)

# dims of PCA to choose. Probably some fine tuning here for dims and number of neighbors, but I think it looks fine enough already
# It is always good to be thinking whether you are overfitting or underfitting data though!
data.merge = FindNeighbors(data.merge, dims = 1:30, reduction = "pca")

Computing nearest neighbor graph
Computing SNN

data.merge = FindClusters(data.merge, resolution = 0.5, cluster.name = "unintegrated_clusters")

Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 5321
Number of edges: 181553

Running Louvain algorithm...

0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

Maximum modularity in 10 random starts: 0.8922
Number of communities: 17
Elapsed time: 0 seconds

data.merge <- RunUMAP(data.merge, dims = 1:30, reduction = "pca", reduction.name = "umap.unintegrated", verbose = FALSE)

Warning: The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
This message will be shown once per session

# Always visualize!
# The plot is showing batch effects, so we definitely need to integrate
DimPlot(data.merge, reduction = "umap.unintegrated", group.by = c("condition", "unintegrated_clusters"), combine = FALSE)

[[1]]

[[2]]

# Now we integrate the two layers using the Harmony algorithm with the PCA we calculated 
data.integrated <- IntegrateLayers(object = data.merge, method = HarmonyIntegration, orig.reduction = "pca", new.reduction = "harmony", normalization.method = "SCT", verbose = FALSE)

Warning: HarmonyMatrix is deprecated and will be removed in the future from the API in the futureWarning: Warning: The parameters do_pca and npcs are deprecated. They will be ignored for this function call and please remove parameters do_pca and npcs and pass to harmony cell_embeddings directly.Warning: Warning: The parameter tau is deprecated. It will be ignored for this function call and please remove parameter tau in future function calls. Advanced users can set value of parameter tau by using parameter .options and function harmony_options().Warning: Warning: The parameter block.size is deprecated. It will be ignored for this function call and please remove parameter block.size in future function calls. Advanced users can set value of parameter block.size by using parameter .options and function harmony_options().Warning: Warning: The parameter max.iter.harmony is replaced with parameter max_iter. It will be ignored for this function call and please use parameter max_iter in future function calls.Warning: Warning: The parameter max.iter.cluster is deprecated. It will be ignored for this function call and please remove parameter max.iter.cluster in future function calls. Advanced users can set value of parameter max.iter.cluster by using parameter .options and function harmony_options().Warning: Warning: The parameter epsilon.cluster is deprecated. It will be ignored for this function call and please remove parameter epsilon.cluster in future function calls. Advanced users can set value of parameter epsilon.cluster by using parameter .options and function harmony_options().Warning: Warning: The parameter epsilon.harmony is deprecated. It will be ignored for this function call and please remove parameter epsilon.harmony in future function calls. If users want to control if harmony would stop early or not, use parameter early_stop. Advanced users can set value of parameter epsilon.harmony by using parameter .options and function harmony_options().

# After integrating we re-find our neighbors and clusters and run UMAP to check our batch effects visually
data.integrated <- FindNeighbors(data.integrated, reduction = "harmony", dims = 1:30)

Computing nearest neighbor graph
Computing SNN

data.integrated <- FindClusters(data.integrated, resolution = 0.5, cluster.name = "harmony_clusters")

Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 5321
Number of edges: 188646

Running Louvain algorithm...

0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

Maximum modularity in 10 random starts: 0.8902
Number of communities: 14
Elapsed time: 0 seconds

data.integrated <- RunUMAP(data.integrated, reduction = "harmony", dims = 1:30, reduction.name = "umap.harmony", verbose = FALSE)

# Much nicer, AR and TOL form well-mixed clusters
DimPlot(data.integrated, reduction = "umap.harmony", group.by = c("condition", "harmony_clusters"), combine = FALSE)

[[1]]

[[2]]

# I would suggest looking into Azimuth instead of SingleR. Better fine labeling
ref <- celldex::ImmGenData()

data.integrated_fine_pred <- SingleR(test = GetAssayData(data.integrated, layer = 'counts'), 
                   ref = ref, 
                   labels = ref$label.fine)

data.integrated_main_pred <- SingleR(test = GetAssayData(data.integrated, layer = 'counts'), 
                   ref = ref, 
                   labels = ref$label.main)
data.integrated$singleR.labels.fine <- data.integrated_fine_pred$labels[match(rownames(data.integrated@meta.data), rownames(data.integrated_fine_pred))]

data.integrated$singleR.labels.main <- data.integrated_main_pred$labels[match(rownames(data.integrated@meta.data), rownames(data.integrated_main_pred))]

DimPlot(data.integrated, reduction = "umap.harmony", group.by = "singleR.labels.main")

DimPlot(data.integrated, reduction = "umap.harmony", group.by = "singleR.labels.fine")

head(sort(table(data.integrated_fine_pred$labels), decreasing=TRUE))

Macrophages (MF.11CLOSER.SALM3)           Monocytes (MO.6C+II-)     T cells (T.8EFF.OT1.D10LIS)    Macrophages (MF.103-11B+24-)         Endothelial cells (BEC)       Macrophages (MF.II+480LO) 
                           1408                             561                             298                             228                             219                             204 

head(sort(table(data.integrated_main_pred$labels), decreasing=TRUE))

      Macrophages           T cells         Monocytes               NKT           B cells Endothelial cells 
             2807               826               352               334               271               220 

# 826 T cells

# Finding more ways to plot our results, checking annotation of cells and pruned cells
plotScoreHeatmap(data.integrated_main_pred)

plotDeltaDistribution(data.integrated_main_pred)

pattern <- "T cells \\(T\\.(8|CD8).*\\)"
# Shorter way to add metadata
data.integrated$is_cd8 <- grepl(pattern, data.integrated@meta.data$singleR.labels.fine)
DimPlot(data.integrated, reduction = "umap.harmony", group.by = "is_cd8")

cd8_subset <- subset(data.integrated, subset = is_cd8 == TRUE)

# Rerun SCTransform on the subset, since the global structure of variance changed 
cd8_subset <- SCTransform(cd8_subset, vars.to.regress = "percent.mt", verbose = FALSE)

Warning: Different cells and/or features from existing assay SCT

# Sub-clustering of CD8 cells
cd8_subset <- RunPCA(cd8_subset, reduction.name = "cd8.pca", verbose = FALSE)

Warning: Key ‘PC_’ taken, using ‘cd8pca_’ instead

cd8_subset <- FindNeighbors(cd8_subset, reduction = "cd8.pca", dims = 1:30)

Computing nearest neighbor graph
Computing SNN

cd8_subset <- FindClusters(cd8_subset, resolution = 0.5, cluster.name = "cd8_clusters")

Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 1061
Number of edges: 43975

Running Louvain algorithm...

0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

Maximum modularity in 10 random starts: 0.6805
Number of communities: 5
Elapsed time: 0 seconds

cd8_subset <- RunUMAP(cd8_subset, dims = 1:30, reduction = "cd8.pca", verbose = FALSE, reduction.name = "cd8.umap")
DimPlot(cd8_subset, reduction = "cd8.umap", group.by = c("condition", "cd8_clusters"), combine = FALSE)

[[1]]

[[2]]

# Finding DE genes
Idents(cd8_subset) <- "cd8_clusters"
cd8_subset <- PrepSCTFindMarkers(cd8_subset)

Found 2 SCT models. Recorrecting SCT counts using minimum median counts: 5574

  |                                                  | 0 % ~calculating  
  |+++++++++++++++++++++++++                         | 50% ~01s          
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=01s  

tcell.cd8.markers <- FindMarkers(cd8_subset, ident.1 = 0, ident.2 = 2)

For a (much!) faster implementation of the Wilcoxon Rank Sum Test,
(default method for FindMarkers) please install the presto package
--------------------------------------------
install.packages('devtools')
devtools::install_github('immunogenomics/presto')
--------------------------------------------
After installation of presto, Seurat will automatically use the more 
efficient implementation (no further action necessary).
This message will be shown once per session

  |                                                  | 0 % ~calculating  
  |+                                                 | 1 % ~15s          
  |+                                                 | 2 % ~14s          
  |++                                                | 3 % ~14s          
  |++                                                | 4 % ~14s          
  |+++                                               | 5 % ~14s          
  |+++                                               | 6 % ~14s          
  |++++                                              | 7 % ~14s          
  |++++                                              | 8 % ~13s          
  |+++++                                             | 9 % ~13s          
  |+++++                                             | 10% ~13s          
  |++++++                                            | 11% ~13s          
  |++++++                                            | 12% ~13s          
  |+++++++                                           | 13% ~13s          
  |+++++++                                           | 14% ~13s          
  |++++++++                                          | 15% ~13s          
  |++++++++                                          | 16% ~13s          
  |+++++++++                                         | 17% ~12s          
  |+++++++++                                         | 18% ~12s          
  |++++++++++                                        | 19% ~12s          
  |++++++++++                                        | 20% ~12s          
  |+++++++++++                                       | 21% ~12s          
  |+++++++++++                                       | 22% ~12s          
  |++++++++++++                                      | 23% ~12s          
  |++++++++++++                                      | 24% ~11s          
  |+++++++++++++                                     | 25% ~11s          
  |+++++++++++++                                     | 26% ~11s          
  |++++++++++++++                                    | 27% ~11s          
  |++++++++++++++                                    | 28% ~11s          
  |+++++++++++++++                                   | 29% ~11s          
  |+++++++++++++++                                   | 30% ~11s          
  |++++++++++++++++                                  | 31% ~10s          
  |++++++++++++++++                                  | 32% ~10s          
  |+++++++++++++++++                                 | 33% ~10s          
  |+++++++++++++++++                                 | 34% ~10s          
  |++++++++++++++++++                                | 35% ~10s          
  |++++++++++++++++++                                | 36% ~10s          
  |+++++++++++++++++++                               | 37% ~09s          
  |+++++++++++++++++++                               | 38% ~09s          
  |++++++++++++++++++++                              | 39% ~09s          
  |++++++++++++++++++++                              | 40% ~09s          
  |+++++++++++++++++++++                             | 41% ~09s          
  |+++++++++++++++++++++                             | 42% ~09s          
  |++++++++++++++++++++++                            | 43% ~09s          
  |++++++++++++++++++++++                            | 44% ~08s          
  |+++++++++++++++++++++++                           | 45% ~08s          
  |+++++++++++++++++++++++                           | 46% ~08s          
  |++++++++++++++++++++++++                          | 47% ~08s          
  |++++++++++++++++++++++++                          | 48% ~08s          
  |+++++++++++++++++++++++++                         | 49% ~08s          
  |+++++++++++++++++++++++++                         | 50% ~07s          
  |++++++++++++++++++++++++++                        | 51% ~07s          
  |++++++++++++++++++++++++++                        | 52% ~07s          
  |+++++++++++++++++++++++++++                       | 53% ~07s          
  |+++++++++++++++++++++++++++                       | 54% ~07s          
  |++++++++++++++++++++++++++++                      | 55% ~07s          
  |++++++++++++++++++++++++++++                      | 56% ~07s          
  |+++++++++++++++++++++++++++++                     | 57% ~06s          
  |+++++++++++++++++++++++++++++                     | 58% ~06s          
  |++++++++++++++++++++++++++++++                    | 59% ~06s          
  |++++++++++++++++++++++++++++++                    | 60% ~06s          
  |+++++++++++++++++++++++++++++++                   | 61% ~06s          
  |+++++++++++++++++++++++++++++++                   | 62% ~06s          
  |++++++++++++++++++++++++++++++++                  | 63% ~06s          
  |++++++++++++++++++++++++++++++++                  | 64% ~05s          
  |+++++++++++++++++++++++++++++++++                 | 65% ~05s          
  |+++++++++++++++++++++++++++++++++                 | 66% ~05s          
  |++++++++++++++++++++++++++++++++++                | 67% ~05s          
  |++++++++++++++++++++++++++++++++++                | 68% ~05s          
  |+++++++++++++++++++++++++++++++++++               | 69% ~05s          
  |+++++++++++++++++++++++++++++++++++               | 70% ~04s          
  |++++++++++++++++++++++++++++++++++++              | 71% ~04s          
  |++++++++++++++++++++++++++++++++++++              | 72% ~04s          
  |+++++++++++++++++++++++++++++++++++++             | 73% ~04s          
  |+++++++++++++++++++++++++++++++++++++             | 74% ~04s          
  |++++++++++++++++++++++++++++++++++++++            | 75% ~04s          
  |++++++++++++++++++++++++++++++++++++++            | 76% ~04s          
  |+++++++++++++++++++++++++++++++++++++++           | 77% ~03s          
  |+++++++++++++++++++++++++++++++++++++++           | 78% ~03s          
  |++++++++++++++++++++++++++++++++++++++++          | 79% ~03s          
  |++++++++++++++++++++++++++++++++++++++++          | 80% ~03s          
  |+++++++++++++++++++++++++++++++++++++++++         | 81% ~03s          
  |+++++++++++++++++++++++++++++++++++++++++         | 82% ~03s          
  |++++++++++++++++++++++++++++++++++++++++++        | 83% ~03s          
  |++++++++++++++++++++++++++++++++++++++++++        | 84% ~02s          
  |+++++++++++++++++++++++++++++++++++++++++++       | 85% ~02s          
  |+++++++++++++++++++++++++++++++++++++++++++       | 86% ~02s          
  |++++++++++++++++++++++++++++++++++++++++++++      | 87% ~02s          
  |++++++++++++++++++++++++++++++++++++++++++++      | 88% ~02s          
  |+++++++++++++++++++++++++++++++++++++++++++++     | 89% ~02s          
  |+++++++++++++++++++++++++++++++++++++++++++++     | 90% ~01s          
  |++++++++++++++++++++++++++++++++++++++++++++++    | 91% ~01s          
  |++++++++++++++++++++++++++++++++++++++++++++++    | 92% ~01s          
  |+++++++++++++++++++++++++++++++++++++++++++++++   | 93% ~01s          
  |+++++++++++++++++++++++++++++++++++++++++++++++   | 94% ~01s          
  |++++++++++++++++++++++++++++++++++++++++++++++++  | 95% ~01s          
  |++++++++++++++++++++++++++++++++++++++++++++++++  | 96% ~01s          
  |+++++++++++++++++++++++++++++++++++++++++++++++++ | 97% ~00s          
  |+++++++++++++++++++++++++++++++++++++++++++++++++ | 98% ~00s          
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 99% ~00s          
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=15s  

tcell.cd8.markers

FeaturePlot(cd8_subset, features = "Gcnt1", reduction = "cd8.umap")

FeaturePlot(cd8_subset, features = "Gzmb", reduction = "cd8.umap")

FeaturePlot(cd8_subset, features = "Spn", reduction = "cd8.umap")

FeaturePlot(cd8_subset, features = "Sell", reduction = "cd8.umap")

FeaturePlot(cd8_subset, features = "Prf1", reduction = "cd8.umap")

FeaturePlot(cd8_subset, features = "Spn", reduction = "cd8.umap")

FeaturePlot(cd8_subset, features = "Cxcr3", reduction = "cd8.umap")

FeaturePlot(cd8_subset, features = "Havcr2", reduction = "cd8.umap")

FeaturePlot(cd8_subset, features = "Ifng", reduction = "cd8.umap")

FeaturePlot(cd8_subset, features = "Icos", reduction = "cd8.umap")

FeaturePlot(cd8_subset, features = "Prdm1", reduction = "cd8.umap")

FeaturePlot(cd8_subset, features = "Tox", reduction = "cd8.umap")

FeaturePlot(cd8_subset, features = "Gzmk", reduction = "cd8.umap")

#FeaturePlot(cd8_subset, features = "Gbma", reduction = "cd8.umap") NOT FOUND
#FeaturePlot(cd8_subset, features = "Gzmh", reduction = "cd8.umap") NOT FOUND

FeaturePlot(cd8_subset, features = "Gzmb", reduction = "cd8.umap")

FeaturePlot(cd8_subset, features = "Lag3", reduction = "cd8.umap")

FeaturePlot(cd8_subset, features = "Tigit", reduction = "cd8.umap")

FeaturePlot(cd8_subset, features = "Irf4", reduction = "cd8.umap")

FeaturePlot(cd8_subset, features = "Tcf7", reduction = "cd8.umap")

FeaturePlot(cd8_subset, features = "Pdcd1", reduction = "cd8.umap")

FeaturePlot(cd8_subset, features = "Tnf", reduction = "cd8.umap")

FeaturePlot(cd8_subset, features = "Tbx21", reduction = "cd8.umap")

#Save the plots as PNGs

library(ggplot2)

List genes (remove genes that were not found)

genes <- c(“Gcnt1”, “Gzmb”, “Spn”, “Sell”, “Prf1”, “Cxcr3”, “Havcr2”, “Ifng”, “Icos”, “Prdm1”, “Tox”, “Gzmk”, “Lag3”, “Tigit”, “Irf4”, “Tcf7”, “Pdcd1”, “Tnf”, “Tbx21”)

Loop through the genes and create/save plots

for (gene in genes) { plot <- FeaturePlot(cd8_subset, features = gene, reduction = “cd8.umap”) ggsave(filename = paste0(“CD8_”, gene, “.png”), plot = plot, width = 6, height = 4, dpi = 300) }

# From the list Dr. Krummey sent, find the differential expression of the genes of interest (done)

# Look back to just NON-CD4 T cells to see if there is a difference in that and solely using/looking at CD8+ T cells